#!/usr/bin/python
# -*- coding:utf-8 -*-

import numpy
from random import random
import math
from itertools import izip

def PCA(X,n):
    '''
    主成分分析。

    nには次元数(>=1)もしくは、累積寄与率の最低値(<1)を指定する。

    入力はこんな感じ。各データを並べる。
    #A=[[-1,-2, 1],
    #   [ 1, 0, 2],
    #   [ 2,-1, 3],
    #   [ 2,-1, 2]]
    計算。Aを2次元に圧縮するならこんな感じ。
    #PCA(A,2)

    返り値は、(寄与率ベクトル、次元圧縮結果のデータ)というタプル。
    '''
    _X=numpy.matrix(X)
    u=sum(_X)/float(len(_X))
    X_u=_X-u#[x-u for x in X]
    ls,vs=numpy.linalg.eigh(X_u.T*X_u)
    ls/=sum(ls)
    cs=vs*_X.T

    if n<1:
        #累積寄与率の最小値指定
        rr=0
        for i in xrange(1,len(ls)+1):
            rr+=ls[-i]
            if rr>n: break
        ls=ls[-i:]
        cs=cs[-i:]

    else:
        #次元数指定
        ls=ls[-n:]
        cs=cs[-n:]

    cs=cs.T

    if isinstance(X,list):
        return ls.tolist(),cs.tolist()
    else:
        return ls,cs



def __difcost(a,b):
    c=a-b
    return numpy.multiply(c,c).sum()

def __NMF(A,initW,initH,featureL,maxiter):
    W=initW
    H=initH

    for i in xrange(maxiter):
        WH=W*H

        if __difcost(WH,A)==0: break
        #print __difcost(WH,A)
        hn=W.T*A
        hd=W.T*WH

        H=numpy.multiply(H,hn)/hd

        wn=A*H.T
        wd=W*H*H.T

        W=numpy.multiply(W,wn)/wd

    return W,H


def NMF(A,featureL=3,maxiter=1000):
    _A=numpy.matrix(A)
    rowL,colL=numpy.shape(_A)

    W=numpy.matrix([[random() for j in xrange(featureL)] for i in xrange(rowL)])
    H=numpy.matrix([[random() for j in xrange(colL)] for i in xrange(featureL)])

    W,H=__NMF(_A,W,H,featureL,maxiter)

    if isinstance(A,list):
        W=W.tolist() 
        H=H.tolist()

    return W,H

def NMF_SVD(A,featureL=3,maxiter=1000):
    '''
    参考:
    山口,桂吾 堀田,政二 宮原,末治: 非負行列因子分解の初期値設定法とその応用   電子情報通信学会論文誌. D-II, 情報・システム, II-パターン処理	09151923	社団法人電子情報通信学会	20040301	87	3	923-928
    '''
    _A=numpy.matrix(A)
    rowL,colL=numpy.shape(_A)

    if featureL>rowL or featureL>colL:
        raise Exception("yi01### feature should be lower than rank(A) ###")

    U,s,V=numpy.linalg.svd(_A)

    W=numpy.abs(U.T[:featureL].T)
    H=numpy.abs(V.T[:featureL])

    W,H=__NMF(_A,W,H,featureL,maxiter)

    if isinstance(A,list):
        W=W.tolist() 
        H=H.tolist()

    return W,H

def LSI(A,n=0.99):
    '''Latent Semantic Indexing'''
    
    _A=numpy.matrix(A)
    u,s,v=numpy.linalg.svd(_A)
    if n<1:
        ssum=sum(s)
        sd=s/ssum
        for i in xrange(len(sd)):
            if sum(sd[:i+1])>n:
                n=i+1
                break
            
    u=numpy.matrix(u)[:,:n]
    S=numpy.matrix(numpy.diag(s[:n]))
    v=numpy.matrix(v)[:n,:]

    _A=u*S*v
    if isinstance(A,list): _A=_A.tolist()
    return _A

def CosSim(x,y):
    return sum([xi*yi for xi,yi in izip(x,y)])/math.sqrt(sum([xi*xi for xi in x])*sum([yi*yi for yi in y]))
        
#if __name__=="__main__":
    #a=numpy.matrix([[1,2,3,4,5,6,7],
                    #[1,2,3,5,67,7,54],
                    #[1,2,3,5,67,7,54],
                    #[1,2,13,5,67,7,46],
                    #[1,2,3,5,67,7,546],
                    #[1,2,3,5,67,7,56],
                    #[1,2,3,51,67,7,546],
                    #[1,2,3,123,4,41,123],
                    #[1,2,3,234,321,3,321]])
    #w,h=NMF(a,7,5000)
    #print (w*h).tolist()
    